{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Cleavage Site and TAP Prediction </h1>\n",
    "\n",
    "\n",
    "This tutorial illustrates the use of Fred2 to predict the steps of the HLA-I antigen processing pathway including proteasomal cleavage and TAP transport. Fred2 offers a long list of prediction methods and was  designed in such a way that extending Fred2 with your favorite method is easy.\n",
    "\n",
    "This tutorial will entail:\n",
    "- Simple cleavage site/fragment prediction from a list of peptide sequences and protein sequences\n",
    "- Simple TAP prediction methods\n",
    "- Consensus prediction of proteasomal cleavage, TAP, and HLA binding to model the complete antigen processing pathway\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from Fred2.Core import Protein, Peptide, Allele\n",
    "from Fred2.Core import generate_peptides_from_proteins\n",
    "from Fred2.IO import read_fasta\n",
    "from Fred2.CleavagePrediction import CleavageSitePredictorFactory, CleavageFragmentPredictorFactory\n",
    "from Fred2.TAPPrediction import TAPPredictorFactory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h2> Chapter 1: Cleavage Prediction </h2>\n",
    "<br/>\n",
    "Fred2 offers a comprehensive list of state-of-the art proteasomal cleavage prediction methods. Usually one distinguishes cleavage methods into cleavage site and cleavage fragment prediction. Cleavage site prediction methods predict each possible cleavage site within a given amino acid sequence, whereas cleavage fragment prediction methods predict the likelihood of a peptide fragment of being a product of proteasomal cleavage. Additionally, the source of the training data are often also quite different. The majority of prediction tools was train on in vitro data of up to three fully analyzed proteins and distinguish between constitutive and immunoproteasomal cleavage. Some prediction methods use natural HLA ligands as training data as they are products of antigen processing pathway and therefore also of cleavage events.\n",
    "<br/><br/>\n",
    "But all methods start with reading in protein sequences. Fred2 offers several ways of defining Proteins. We can either directly initialize a `Fred2.Core.Protein` object by specifying a amino acid sequence and optionally a progenitor gene and transcript id, as well as the progenitor `Fred2.Core.Transcript` object, or we can directly read in proteins from FASTA files by using `Fred2.IO.read_fasta`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PROTEIN: AAAAAAAAAAA (aa-seq)\n",
       "\t\t  Dummy (orig transcript)\n",
       "\t\t VARIANTS:"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "protein = Protein(\"AAAAAAAAAAA\", gene_id=\"Dummy\", transcript_id=\"Dummy\")\n",
    "proteins = read_fasta(\"./data/proteins.fasta\", id_position=3, in_type=Protein)\n",
    "protein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have a protein sequence to work with, we can specify the cleavage site prediction method of our choice. Fred2 offers one entry point for each type of prediction methods via so called factories. For cleavage site prediction it is `CleavageSitePredictorFactory`. To get an overview which prediction methods are currently implemented, we can use `CleavageSitePredictorFactory` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proteasmm_i 1.0\n",
      "netchop 3.1\n",
      "proteasmm_c 1.0\n",
      "pcm 1.0\n"
     ]
    }
   ],
   "source": [
    "for name,version in CleavageSitePredictorFactory.available_methods().iteritems():\n",
    "    print name, \",\".join(version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets select `PCM` for example and make predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Seq</th>\n",
       "      <th>pcm</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th>Pos</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">Protein_0</th>\n",
       "      <th>0</th>\n",
       "      <td>M</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>K</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>M</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>R</td>\n",
       "      <td>-1.77936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>F</td>\n",
       "      <td>-0.87673</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Seq      pcm\n",
       "ID        Pos             \n",
       "Protein_0 0     M  0.00000\n",
       "          1     K  0.00000\n",
       "          2     M  0.00000\n",
       "          3     R -1.77936\n",
       "          4     F -0.87673"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcm = CleavageSitePredictorFactory(\"PCM\")\n",
    "site_result = pcm.predict(proteins)\n",
    "site_result.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To specify a particular version of a prediction method, we can use the flag `version=\"\"` when calling the PredictorFactories. If we do not specify any version, Fred2 will initialize the most recent version that is supported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Seq</th>\n",
       "      <th>pcm</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th>Pos</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">Protein_0</th>\n",
       "      <th>0</th>\n",
       "      <td>M</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>K</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>M</td>\n",
       "      <td>0.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>R</td>\n",
       "      <td>-1.77936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>F</td>\n",
       "      <td>-0.87673</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Seq      pcm\n",
       "ID        Pos             \n",
       "Protein_0 0     M  0.00000\n",
       "          1     K  0.00000\n",
       "          2     M  0.00000\n",
       "          3     R -1.77936\n",
       "          4     F -0.87673"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcm = CleavageSitePredictorFactory(\"PCM\", version=\"1.0\")\n",
    "site_result = pcm.predict(proteins)\n",
    "site_result.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "External tools like `NetChop` offer two additional flags when calling `.predict()`, `command=\"/path/to/binary\"` and `options=\"command options\"`. `command=\"\"` specifies the path to an alternative binary that should be used instead of the one which is globally registered. With options you can specify additional commands that will directly be passed to the command line call without any sanity checks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For CleavageFragment prediction we first have to generate peptides from protein sequences with `Fred2.Core.generate_peptides_from_proteins`. Fred2 currently supports only one CleavageFragment prediction methods proposed by Ginodi et al. <a href=\"http://bioinformatics.oxfordjournals.org/content/24/4/477.long\">(Ginodi, et al.(2008) Bioinformatics 24(4))</a>, which supports only 11mers (9mer epitopes and two flaking amino acids)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ginodi</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(A, A, V, D, P, A, D, R, C)</th>\n",
       "      <td>-18.891867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, D, R, C, K, E, V, Q, Q)</th>\n",
       "      <td>-14.992024</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, V, D, P, A, D, R, C, K)</th>\n",
       "      <td>-14.755080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, G, K, A, A, V, D, P, A)</th>\n",
       "      <td>-17.551360</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, K, E, V, Q, Q, I, R, D)</th>\n",
       "      <td>-17.682320</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                ginodi\n",
       "(A, A, V, D, P, A, D, R, C) -18.891867\n",
       "(A, D, R, C, K, E, V, Q, Q) -14.992024\n",
       "(A, V, D, P, A, D, R, C, K) -14.755080\n",
       "(C, G, K, A, A, V, D, P, A) -17.551360\n",
       "(C, K, E, V, Q, Q, I, R, D) -17.682320"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pep = generate_peptides_from_proteins(proteins, 11)\n",
    "CleavageFragmentPredictorFactory(\"Ginodi\").predict(pep).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result object is based on pandas' `DataFrame`, thus all possibilities of manipulating the results pandas offers are possible, including rudimentary plotting capabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEPCAYAAABiCi5wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XvUJVV95vHvIyAQwQiECSBNQAWWIjLAiLqighCg0ay0\nDEOEcQwXL6gLVCCjAit2G0SDGdEl6DiJtIkmXESJ0kFGGKC9gaLSNE2n5drIRQcb0QAKAza/+aPq\nherz1qmzT1Wd+/NZ66z3XH5n7127qvZbZ1ftXYoIzMxscjxr1AUwM7P+uOE2M5swbrjNzCaMG24z\nswnjhtvMbMK44TYzmzCNG25JSyU9IGlVGwUyM7NqbRxxfwFY2EI6ZmaWoHHDHRHfAX7VQlnMzCyB\n+7jNzCaMG24zswmz8aAzkOTJUMzMaogIdfug8QPYGVjV5bMoeW9JYrpjGzfOZfOyTk6dAFH2mMZl\nnaa4YeRZtR20cTnghcB1wG6S7pV0XNM0zWZLAIt5pt02q9a4qyQijm6jIGZmlmbgfdxdLJ+CuFHk\nOaq4UeTZSlznORZJi6Gi7zA935SYPuIOSEzO639M4kaR59OU96UMjKTosZOYbaDshHbdbShLqzM5\n1U6vbeNePhudqrZzVEfcZj0UGzO3YWZFbrj71O3yRh8hmQ3PrO+HbVxVslDSTyTdLukDbRRq/HVe\nvWVmwze7+2GjhlvSRsB5ZJNMvQQ4WtKL2yiYmZmVa3rEvR9wR0TcHRFPAhcBi5oXy8zMumnax/18\n4N7C6/uAVzRM02zoUvtMZ71v1cZD0yPupI4lSUvyx3JJ0fFY0iW+35hucZ0xUbbzpcbl0R2P0mVu\nLd+2l2Hc4/Jo2qjf+WkNp31ts3zjvr6GEVen3ka8DKntU1nc8rLYDb7X5DpuSa8kG2O/MH99GvBU\nRJxdXNBRHo1kldr7OtnUuFHk2/YyjHtcirbXV9tGsazTGzce+2E/cW2oajubHnH/CNhV0s6Sng28\nCbisYZpmZlahUR93RPxO0onAN4GNgPMjYk0rJTMzG0uj/1HXxiRTVwBXtFAWM7MBaK+hHZeuOI+c\nNLOpNS4Nbdt86zIzswnjI24zG6KpPAAeOjfcZjYU09ptMQpN5yo5UtJqSesl7dNWoabH8Ad+jD/X\niVlTTY+4VwGHA/+rhbJMFR9dzNd+ncxSFc/SslovTa/j/gmANC0b1bQsx/SbpX+Ms7SslsZ93Dnv\nHGY2KXo23JKuArYr+ej0iFiWkknHhCnLI2J5UunMzGaEpANIvGt0KzcLlnQtcGpE3Fjy2URMMjUq\nszTJ1CzpNqPdIOtkWtZ/m5NM1c9zOPlWqWo72+wqmdmd1KzTLP/TssFrejng4ZLuBV4JXC7Jc5aY\nmQ1YK10llRm4q6SSu0qsLdOy/t1Vkuc+wPm4zcxsyNxwm5lNGDfcZmYTpunJyb+VtEbSSkmXSvr9\ntgpmZmblmh5xXwnsERF7AbcBpzUvkpmND08KNo4aNdwRcVVEPJW//AGwY/MiDYI3PrN+RYTKHqMu\nl7U7AOd44MIW02uFNzQzmzatzFUi6QzgiYi4oOXymZlZh54Nd0QcXPW5pGOB1wMHVcQsKbz0JFNm\nZh2GNsmUpIXAJ4D9I+LBLjEjHTk57jxy0toyPiMd207PIyc7Nb2q5FxgC+AqSSskfbZhemZm1kPT\nO+Ds2lZBzGyY2j6Q9I+rYfIdcMxmTNs//90lNnxuuK0H75Nm48YNt3XlIymz8eSG28wMmKRfl7Wv\nKpF0Zj651E2Srpa0oM2CmZkNy6QN7699HbekLSPikfz5ScBeEfG2kjhfx11hnK/jbnMZbPCmZT2M\n4jrucTSQ67jnGu3cFkDpABwzM2tXoz5uSWcBbwF+S3bDYDMzG7DKrpKUCabyuA8Cu0fEcSVpBPDh\nwlueq6TAXSXWlmlZD7PaVVIyV8nibsvdyl3eJe0EfCMiXlrymfu4K7jhtrZMy3qY1Ya700D6uCUV\nh7svAlbUTcvMzNI16eP+mKTdgfXAncC72imSmdU3cwemM6mVrpLKDNxVUsldJWYbcldJZpDTupqZ\n2ZC54TYzmzBuuM3MJownmTKzMTRzXdp9aXzELelUSU9J2rqNApnZbJukyZ5GpVHDnc8IeDDw03aK\nY2ZmvTQ94j4HeH8bBTEzszRNRk4uAu6LiJtbLI+ZmfVQeXKyYpKpM4DTgEOK4RXpLCm89CRTZmYd\nSiaZ6h5bZ+SkpJcCV5NN5wqwI3A/sF9E/KIj1iMnK3jkpJmVqWo725odcC2wb0Q81E/m5obbzMoN\nY8j7YCc8MTOzp7UyACciXtBGOmZm1puHvJuZTRg33GZmE8YNt5nZhGkyAGeJpPskrcgfC9ssmJmZ\nlWtycjKAcyLinLYKY2azIbuE1Ob0e+ls06tKfJ2umdXi6/wzdf6JNe3jPknSSknnS3pew7TMzCxB\n5cjJHnOVfB9Yl78+E9g+It5akoZHTlZof+TkfB45aePG7cIzutVFVR1VdpVExMGJGX8eWFbx+ZLC\nS08yNSDeEcwm18Anmcoz2T4ifp4/Pxl4eUT815I4/2et0OYRd5t59p/efF7v1o3bhWe0fsTdw9mS\n/iNZC7AWOKFBWjbBvAOaDVcrswNWZuD/rJWm4YjbrF9l7cIwLhEcx2182EfcZmYtG2TbPXZtdm0e\n8j4W1PEws1GQdLekD0paLekhSUslbZp/tkjSTZL+XdIdkg7J318u6UxJ35P0iKTLJP2BpH/OY2+Q\n9EetltNdJePPXSU2bbp3lQz2iLvXNi7pbuBh4DCyO3wtA67N/14JHBERV0vaAdgyIm6VtBzYATgU\n+CVwPfBssvN+3wKWAusj4vguebqrxMysgQDOi4j7ASSdBZwLbAucHxFXA0TEzzq+84WIWJt/5wrg\nxRFxTf76ErKxLq1p1FUi6SRJayTdIunstgplZjZC9xae30N2NL0AuKviOw8Unj8O/KLj9RatlY4G\nR9ySXgf8GfCyiHhS0rbtFcvMbGR26nj+M7LG/IWJ3x/41TFNjrjfBXwsIp4EiIh1PeLNzMadgHdL\ner6krcmm97gIOB84TtKBkp6Vf757x/fKng9Ek4Z7V+C1kr6fn1X9T20VyobBV7LYOOrcLtt8JAng\nArITkXcCtwMfiYgfAscBnwR+DSxnwyPz6HjeedTd6lF4k0mmzgKuiYj3Sno5cHHZTYPzi+o/XHjL\nc5X0yVeB2LQZ16vNJK0F3jp3YnFIeUZEqGSuksXd6qjJXCVXAH8TEd/KX98BvCIifllWqFqZGOCG\n26bPuLYLo2y4U9+HZl0lXwMOzDPYDXh2Z6NtZmbta3Id91JgqaRVwBPAX7RTJDOz0YiIXUZdhhQe\nOTkB3FVi08btwjOG3VViZmYj4IbbzGzCeK6SieFflTZdhjH/9rRywz0B3Bdo08bbdDO1u0okXSRp\nRf5YK2lFH989YNLjxrlsbceNc9lGFTfOZWs7bpzLNqq4UZVtTu2GOyKOioi9I2Jv4Kv5I9UBUxA3\nijxHFTeKPMc9bhR5jipuFHmOe9wo8nxa464SSQL+HHhd07TMzKy3Nq4qeQ3wQETc2UJaZmbWQ91J\npk6PiGV5zP8EbouIT3ZJw2eOzcxqaH2SKQBJGwP3Aft03MrHzMwGpGlXyZ8Aa9xom5kNT9OG+03A\nhW0UxMzM0gx8kikzM2vX0OYqkfSUpC8VXm8saZ2kZR1xjxaev17SrZIW1Ekr/2x9PkjoZkmXSiq9\n23Ihbu6xU5O4PPYMSbdIWpnH7lc3rpDvKklflrR5l7S2ywdH3SHpR5Iul7RrR8wfSrpA0p15zHWS\n3liSVs91Ufj8jfl62b3s89Q4SddIOqTjvfdJ+mzHe3P1cYukmySdkl+a2pleUr0llq3fPCu3uTw2\nZX0lrYeSbfP9JTFJ9Zu/v6Okr0u6LS/fpyRtUjeubHkqliF1f+21LlKX4dGO18dKOrci35TtqXK/\n7mf/micihvIAHgFuBDbLXx8GrAAu64zL/x5Edr+3XeqmVUwvf/4PwKndype6HIlxrwKuAzbJX28N\nbN8grrgc/wScXBIj4HrgHYX3Xga8ukfMTsCJ3fKsWheF2IuBy4AlPeqlMg54O7C0473ri8tQUh/b\nAleVpZlSb32UrU6eVdtcz/XVz3pI2Tb7qF8BNwDH5K+fBXwe+HiduNRy9lF3PddFP2XrLBNwDHBu\n3e2JhP26n/1rXvqpgU0fZI3tR4Aj8tdfBN4PLCuJey3ZjTp3a5JWSUWfAHy27kbfZ9zhlPwjaRBX\nXI53Ap8piTkQ+FaPdA4iu+9n6jqrXBd53BbA3WT/ANY0ics38AeAjfPXOwM/7bUegF2AB+vUWx9l\nq5Nn1TbXc331uR5SGu7U+j2os2zAlsCD5AdM/cSllrPu/lq2LvopW0l6x9K74a7annru16nrtewx\n7GldLwaOkrQpsCfwg5KYzYB/ARZFxG0N03qapI2AQ4BbuoRsXviJWTV8PzXuSmBB/vPnM5Je2zBu\nbjk2JvuFcXPJxy8Fflz1fWAPsl8rKVLXxSLgf0fEPcA6SfvUjYuIh8iOkl6fv3UU2bquFBFrgY0k\nbVv2eY9662cZ+smz1zaXsr4gfT0Ut80Vko4sKXNq/e7RWbaIeAS4B3hRjbi+JNTdBrqsi37KtkHd\nkd3gvOsJwITtKWW/Tl2v8wy14Y6IVWT/4Y8GLu8S9gTwPeBtLaQF+QoBfg4sAD7XJe6xyOdeiYgj\nKtJLiouI3wD7Au8A1gEXSzqmblxhOX5IdmR4flm2FeUujZF0Xt5HeENJbNK6IFsHl+TPL8lfN4m7\nkKxBgeZXLqXUWz9l6yfPXttc6pUBqeuhuG3uHRGXdIlLqd+qshWnykiNS5Vadyn6KdsGdQd8iPK5\nlJO2p8T9OnW9lmYwlAfP9Of8FdlPlT3IJlYp6yrZjKx/6LQmaXXEbg58Gzg85adXr+WosfxHkNYl\nUhqXki9pXSUH0tFVAmwDrC3LM2FdbA38hmwjXkt2NFP20zspLo/dguzn/N7ArSnrAXgBPbotKuok\ndRn6yjNhm+unq6RyPfS5DafUb1k3w3PJGqHn9BuXWs66+2vZuuinbCXpHUuPrpJ+HmX7dep6LXuM\n4g44S8lOIqzuFhARjwNvAN4s6fgmaRXSfAx4D3BW2dnntknarePqgL3JGoZacSki4hpgU0lvL6T/\nMkmv7ojZTNI7C199TkWavdbFfwG+GBE7R8QuEbETsFbSa2rGERGPAtcCXwAu6LXc+c/jzwHzrgJI\nlFy2fvLstc2lrK9CbOo+0VNK/UbE1cDvSXpLXq6NgE8AF0R2NNlXXI0yJu+v3dbFoMqWInW/rr1e\n6/z3qPkf5+GS9/Zn/n+hhwvPdwTuAv60TlplsWRXDbwppXypy9Elbh+yn0GrgZXAV4CtG8Sl5rs9\nWZ/lHWT9g8uAF3bEbEf28/gusnMD1wBHVuVZsS6uAQ7peO8kOk7apMYVPlsErKf7CerfkV1JdAtw\nE3AK+biEfuutj2WolWe3ba6P9dVzPXSUb+7x0YplrqzfQl5fB24DfpWXc5O6cSnrpI/9NXVdpC5D\nZ77HAJ+usz3lcT3369T1WvbwABwz60nSq4C/J/sHv6Zp3CiMc9n65YbbzGzC+C7vZmYTpnHDLWmB\npGslrc6Hd76njYKZmVm5xl0lkrYDtouIm/J5BX4MvHHS+5DMzMZV4yPuiPi/EXFT/vxRYA2wQ9N0\nzcysXKt93JJ2JrtesXL4uZmZ1df4Lu9z8m6SrwDvzY+85973ZStmZjVEl3tO9rzQO/Fi802AbwLv\nK/ksSt5bkpBmKzHOb/rLVBXHhrO53UzFtK7junzjWKZpz28cylTWds492riqRGQTrfxbRHyqaXpm\nA3QP8CJJW0n6Wj7B/fWS9gSQtD9wQj5D3I2quAGC2Si10cf9x8B/A15XmBZxYQvpmrUmn4bzRWRH\n3X8N/Dgi9gJOJ5vPHeBU4PLIZod7NfDYKMpq1kvjPu6I+C79/wNYPsQY5zf8mHHKb24aTsjmtVhK\ndvL8PwNExLWStpG0JdncEsdIOgm4NCLub1iuYcZMZX6d58gkLYaKvt8prIMyAx/yLikqKtlsoCQ9\nEhFbdrx3I9ndk9bmr+8BXhIRj0rag2y2tncDh0bErUMvtD0ta7g72yhVNdxTo6rt9JB3m0XfAd4M\nIOkAYF3eaL8wIlZHxMfJJsqvvPGx2ai0cXJyqaQHJK1qo0BmLSv7SbkE2FfSSuCjZFN4ArxX2d27\nV5LdneSK4RRx+kiKzseoyzRN2hjy/hrgUbKJ6Pcs+dxdJWYzZn4XR73uDXeVDKirJCK+QzZBuZmZ\nDYH7uM3MJkxrQ97NzCZdt774ceuaGUrDLWlJ4eXyiFg+jHzNzPo3v099GPIrnA5Iim3jOu58VsBl\nPjlpZjC5JyfH6WToQE9OSroQuA7YTdK9ko5rmqaZDV/ZJXy+jG88eeSkmQHtHm36iLudsnTL1ycn\nzcbUpJwos+Fzw2021to5UVb2T8D/ACZXG33cCyX9RNLtkj7QRqHMbBCi8LBJ1qjhlrQRcB6wEHgJ\ncLSkF7dRMDMzK9f0iHs/4I6IuDsingQuAhY1L5aZmXXTtI/7+cC9hdf3Aa9omGYrfGJnsnn9mXXX\ntOEe886y0YyAsrZ4/ZmVadpVcj+woPB6AdlR9wYkLckfyzsu7l/SJbbfmNK4bEcvPkoiEgYcpAxM\nSB28MCn5jbpMeSRV62/a67zNOpifVjfDyy89rcHtx73L3Wi7q9uWLS+L3eB7TQbgKLsB663AQcDP\ngBuAoyNiTXEhx/nnbVbh1QMF5sfMj0uJmZz8Rl+mFNNe5ylSy9SWNvMbbFqTX+dVbWejrpKI+J2k\nE4FvAhsB5xcbbTMza18bd3m/At/iyUZmbH/MmQ2MR04C3vkn0zh3wdm0G+2mN/MNd/rO7zbCzMbj\ngGHmG+4U47CizMzmNB3yfqSk1ZLWS9qnrUKZmVl3Ta/jXgUcDny7hbKYmTWQcp36dGh6OeBPAKTp\nrygzG1+z1p3pPm4bkZnaz8xa1bPhlnQVsF3JR6dHxLKUTOS7vFvBbBwdzcAiWqs0gru8XwucGhE3\nlnw21kPe2zLtw6/bLNOwlc0jMerlm7Th14PIb/RlH/22WWVgQ94782kxrQnlKhhHvlbfpk3TywEP\nl3Qv8ErgckkzO/Q9IlT2GHW5LI3Xn02SVrpKKjOYka6SVO4qmVzuKinvdoJ65y3cVVJtWF0lZjbl\nxrmhmyWN7/JuZmbD1bSP+28lrZG0UtKlkn6/rYKZmVm5pkfcVwJ7RMRewG3Aac2LZGZmVRo13BFx\nVUQ8lb/8AbBj8yKZzYrZmVvD2tXmycnjgQtbTM9sDLXTyPoknzXRypB3SWcAT0TEBV3SWFJ46SHv\nNpHc2NogDXXIu6RjgbcDB0XE4yWf+zruAl/HbW2a5Dr3ddzVBnYdt6SFwH8H9i9rtM3Mqo1tuznW\nGh1xS7odeDbwUP7W9RHx7o4YH3EX+Ijb2uQ6T+cj7lxE7Nrk+2ZmwzW27XRfPOTdzGbCOB9d98tD\n3s3MJkzthlvSmflQ95skXS1pQZsFMzOzcrVPTkraMiIeyZ+fBOwVEW8rifPJyQKfnLQ2tTnNqo2X\ngZycnGu0c1sAD9ZNy8zqcQM9m5pex30W8Bbgt2R3wTEzswGr7CpJvcO7pA8Cu0fEcSVpuKukYPq7\nSubz+jfrX+2ukog4ODGPC4BvVBRgSeGl5yqZUm6gzeobylwlknaNiNvz5ycB+0XEW0rifMRdMM1H\n3GbWnkGNnPyYpN2B9cCdwLsapGVmZol8l/ch8xG3maWoajs9ctLMbMK44TYzmzBuuM3MJkzjhlvS\nqZKekrR1GwUyM7NqjRrufGKpg4GftlMcMzPrpekR9znA+9soiJmZpWkyresi4L6IuLnF8thYU8fD\nzEahcgBOxVwlZwCnAYcUwyvSWVJ46SHvE8jXa5sN1sCHvEt6KXA12ayAADsC95MNe/9FR6wH4BRM\n6gAcMxuuqrazlZGTktYC+0bEQyWfueEucMNtZimGMXJysOPmzczsaa3c5T0iXtBGOmZm1lsrDbeZ\nWT+63XRjVvXbNemG28xGwudRMnX+iTW5jnuJpPskrcgfC+umZWZm6ZoccQdwTkSc01ZhrMgHI2ZW\nrmlXiVuXAfBPSDOr0vRywJMkrZR0vqTntVIiMzOrVDkAp8eQ9+8D6/LXZwLbR8RbS9II4MOFt2Z6\nyHvZiYg6R9gegGOTzAPznjFXFyVD3hcPeuTkzsCyiNizW6EaZ2IbcMNtk6ysXRjGJYLjuF13ayMH\nMnJS0vaFl4cDq+qmZWaWiQE+pkeTPu6zJd0saSWwP3ByS2WyZJ5m1axNku6W9EFJqyU9JGmppE3z\nzxZJuknSv0u6Q9Ih+fvLJZ0p6XuSHpF0maQ/kPTPeewNkv6o1XK20VVSmYG7SkYqpU/dXSU2bN27\nSgbZHvXeriXdDTwMHEY2++ky4Nr875XAERFxtaQdgC0j4lZJy4EdgEOBXwLXA88GTgC+BSwF1kfE\n8V3y7LurxCMnp5wbYLO+BHBeRNwPIOks4FxgW+D8iLgaICJ+1vGdL0TE2vw7VwAvjohr8teXkF3A\n0Rrf5d3MbEP3Fp7fQ3Y0vQC4q+I7DxSePw78ouP1Fq2VjuY3Cz5J0hpJt0g6u61CmZmN0E4dz39G\n1pi/MPH7Az8TWrurRNLrgD8DXhYRT0ratr1imZmNhIB3S/pX4DGyMSsXAV8HrszfXw5sD2wREbcW\nvkfJ84FocsT9LuBjEfEkQESs6xFvZtZD55VSbT6SBHAB2YnIO4HbgY9ExA+B44BPAr8ma7x36vhe\n8XnnUXerR+G1ryqRtILsv9BCsj6cv4yIH5XE+aqSMeerSmzYxrVdyG/D+Na5E4tDyrPdq0p6DHnf\nGNgqIl4p6eXAlwHfCcfMbMAqG+6IOLjbZ5LeBVyax/1Q0lOStomIX5bELim8nOm5SsbX2B38mM2U\nkrlKusc26Co5AdghIhZL2g34PxGxU0ncWP4kMrPRcbvwjGEPwFkKLJW0CngC+IsGaZmZWSIPeTez\noXO78Iyhzg5oZmaj4blKzGwkhjH/9rRqMnLyImD3/OXzgF9HxN6tlMrMppq7SZqp3VUSEUdFxN55\nY/3V/JEkv+xlKDHOb/rLNO35jWOZpj2/cSxTUeM+bkkC/hy4sI+vHTDEGOc3/BjnN/wY5zf8mFHk\nB7RzcvI1wAMRcWcLaZmZWQ91h7yfHhHL8udHk03KYmZmQ9DoOm5JGwP3Aft03BGiGOMzx2ZmNQzq\n1mV/Aqzp1mhXZWxmZvU07eN+E/2dlDQzs4YGPuTdzMzaNfQh75Ie7fH5ekkrJK2S9GVJm3d8fo2k\nQzree5+kz9bJLyWuUKa5x/srYm6WdKmkeTcHVeY7khYW3jsyvyv0vHJIer2kWyUtKEnrDyVdIOlO\nST+SdJ2kN/Zb9j7y21HS1yXdJukOSZ+StElJXM/6bmsd95nXLZJuknRKfglraTo96mCbQj3+XNJ9\n+fMb5+qij7Q610vZzJq99pXkfSGhznvWU1mZJB0r6dw65VI2FfT/KLz+S0mLOz7/UuH1xpLWSVpG\nB0ln5GVfmS/HfhV1ULUft1LnqdtBar7zRMRQH8AjqZ8D/wSc3PH524GlHe9dD7y6Tn4pcSlpdJT7\nH4BTu8TtAfwbsCnZnZ9vA3bpTAc4iOy2SbuUpKF8md9ReG8n4MS69Z2Q3w3AMfnrZwGfBz7eQl3V\nXsc18toWuApY0m8dlKS7GDilTn3WKXuXz5P3hYQ671lPZWUCjgHOrbn+Hie7Pdg2+etTgcXFvIAb\ngc3y14cBK4DLOtJ5FXAdsEn+emtg+1HWeb/bVErZio9xn2Tqu8CLOt77KvAGZVe0IGlnsnnBvzvc\nolW6ni53hI6I1cAy4APAh4B/jIi1xRhJrwX+DnhD52e5A4H/FxF/V0j3nog4r05hE/N7LCL+Mc/r\nKeBk4HhJm9XJs2Bo6ziy+6K+Azix87OEOihTeuK9Zlp11K2nsjp/WlU9lSirg9RyPUlWTydXpP8N\n4A3586PJzql15rkd8GA8c//bhyLi5wllryO5zge5HYxtw51XzGHAzcX3I+IhsqO/1+dvHQVcPODi\nbN7xE+vIboGSNgIOAW6pSO/DwJuBQ4GPd3y2GfAvwKKIuK3L9/cgOxJpo+yp+f24+EZEPALcQ0UD\n0Mso1nG+A20kadvC2yl1kCo1reJ6SZ4uoqhOPXWr85K0y+qps9wryLbl6PhuP+X6LPBmSc/t8vnF\nwFGSNgX2BH5QEnMlsCDvjvhM3mCWSd6Pu+lj2drcpuYZx4Z783yD+CFwN3B+ScyFZBUGw7my5bHI\n52XJH5eUxMyV++fAAuBz3RKLiN8CFwFfmjtKKHgC+B7wtorybLCjSDov75e8oUbZ+86vQ51LSsdt\nHafUQdtpFdfLEQ3yS62nlDpPscH2RParseyoO6lc+QHAF4H3dPl8FbAz2dH25V1ifgPsS/YrYR1w\nsaRjepW9y36cImXZ2tym5uunX6WNB330cVfEbAE8AOwN3NpG31FVXGKZ5vq0Nge+DRzeI34xJf3g\nZP16m5H12Z3W5bsHkt27s/jeNsDaOvWdkN9BwLc63nsu2U7ynLp11XQd18mL7IbWD/ZbB22uv2HX\nU+p20KueusQdS0cfd7/rD9gKWEv2T2Bxyed/BTxI9svvAGBZj2U5go5+8FFsm/1sUyn5Fh/jeMTd\nU0Q8ClwLfIExG24fEY+RHT2cVXZWPjGNx8n69d4s6fiSz68BNpP0zsLbz6mTV2J+VwO/J+kt8HR3\n0CeACyI72mndINZx/rP/c8C5nZ/1qoN+tJlWQl5DradBlCsifgV8GXgr5b/ulpKdKF3dpby7Sdq1\n8NbeZL8qBiJ12Qa6HfTTyrfxAB5u8nkhbhGwHtitpfS6xgG/IzubPff4aK/vA5cBb6pIs9tVCQ8X\nnu8I3AX8aUncdmQ/0e4i6/e7Bjiy37L3kd+OwNfJroL5FVm/3iYlceuBewuP9w1qHSfmNbf8twA3\nAaeQj19ujPN+AAABxElEQVTotw5aXn896yBl+VL3hYT9rmc9ddnOjwE+XXP9FevqPwC/AT5UVWZg\nf+ZfVbIPWbfEamAl8BVg6373hTbrvN9tKnWfmHt4AI71TdKrgL8n+0exZtTlMZs1brjNzCbMRPZx\nm5nNMjfcZmYTxg23mdmEccNtZjZh3HCbmU2YpnfAMRtrktaTzcuxMbCGbIbDx0ZbKrNmfMRt0+63\nkc1LsSfZ/BHv7PUFs3HnhttmyXeBF0naStLX8kn3r5e0J4Ck/Qszx92okpthmI0DN9w2E/LpTBeS\ndZv8NfDjiNgLOJ1sdjrIJvJ/d2Sz3r0acJeKjSU33DbtitOZ/pRswqI/Br4EEBHXAttI2pJsvotP\nSjoJ2Coi1o+ozGaVfHLSpt1j+RH00/JJGztnboyIOFvSv5LN6PY9SYdGxK1DKqdZMh9x2yz6Dtkd\niJB0ALAuIh6V9MKIWB0RHyc7Qt99hGU068qTTNlUk/RwRDy3472tyLpMXkA2leg7IuIWSZ8GXgc8\nRTa96bEx/w5FZiPnhtvMbMK4q8TMbMK44TYzmzBuuM3MJowbbjOzCeOG28xswrjhNjObMG64zcwm\njBtuM7MJ8/8BEU4Gbkliz/wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11470db90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "f, a = plt.subplots(len(site_result.index.levels[0]),1)\n",
    "\n",
    "for i,r in enumerate(site_result.index.levels[0]):\n",
    "    site_result.xs(r).plot(kind='bar', ax=a[i]).set_xticklabels(site_result.loc[(r,slice(None)),\"Seq\"], rotation=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also can combine several prediction results of the same type via `CleavageSitePredictionResults.merge_results` (Note this function returns a merge Result DataFrame):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Seq</th>\n",
       "      <th>pcm</th>\n",
       "      <th>proteasmm_c</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th>Pos</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"7\" valign=\"top\">Protein_0</th>\n",
       "      <th>0</th>\n",
       "      <td>M</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>K</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>M</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>R</td>\n",
       "      <td>-1.779360</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>F</td>\n",
       "      <td>-0.876730</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>F</td>\n",
       "      <td>0.379716</td>\n",
       "      <td>0.967189</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>S</td>\n",
       "      <td>-0.682109</td>\n",
       "      <td>0.848493</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Seq       pcm  proteasmm_c\n",
       "ID        Pos                           \n",
       "Protein_0 0     M  0.000000     0.000000\n",
       "          1     K  0.000000     0.000000\n",
       "          2     M  0.000000     0.000000\n",
       "          3     R -1.779360     0.000000\n",
       "          4     F -0.876730     0.000000\n",
       "          5     F  0.379716     0.967189\n",
       "          6     S -0.682109     0.848493"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy\n",
    "\n",
    "site_result2 = CleavageSitePredictorFactory(\"proteasmm_c\").predict(proteins)\n",
    "merged_result = site_result.merge_results([site_result2])\n",
    "merged_result.head(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can also filter the results based on multiple expressions with `CleavageSitePredictionResults.filter_result`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Seq</th>\n",
       "      <th>pcm</th>\n",
       "      <th>proteasmm_c</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th>Pos</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"7\" valign=\"top\">Protein_0</th>\n",
       "      <th>5</th>\n",
       "      <td>F</td>\n",
       "      <td>0.379716</td>\n",
       "      <td>0.967189</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>V</td>\n",
       "      <td>0.563343</td>\n",
       "      <td>0.806216</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>P</td>\n",
       "      <td>0.662944</td>\n",
       "      <td>1.011543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>A</td>\n",
       "      <td>0.115836</td>\n",
       "      <td>0.740608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>D</td>\n",
       "      <td>0.124760</td>\n",
       "      <td>1.060129</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>C</td>\n",
       "      <td>0.853636</td>\n",
       "      <td>0.963365</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>D</td>\n",
       "      <td>0.445906</td>\n",
       "      <td>1.130923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"8\" valign=\"top\">Protein_1</th>\n",
       "      <th>3</th>\n",
       "      <td>I</td>\n",
       "      <td>0.159584</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>V</td>\n",
       "      <td>0.168248</td>\n",
       "      <td>1.003839</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>L</td>\n",
       "      <td>0.589015</td>\n",
       "      <td>1.141511</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>P</td>\n",
       "      <td>0.523189</td>\n",
       "      <td>0.928291</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>D</td>\n",
       "      <td>0.804995</td>\n",
       "      <td>1.003516</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>V</td>\n",
       "      <td>1.404717</td>\n",
       "      <td>0.989936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>M</td>\n",
       "      <td>1.662398</td>\n",
       "      <td>1.064119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>L</td>\n",
       "      <td>1.071719</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Seq       pcm  proteasmm_c\n",
       "ID        Pos                           \n",
       "Protein_0 5     F  0.379716     0.967189\n",
       "          14    V  0.563343     0.806216\n",
       "          16    P  0.662944     1.011543\n",
       "          17    A  0.115836     0.740608\n",
       "          18    D  0.124760     1.060129\n",
       "          20    C  0.853636     0.963365\n",
       "          28    D  0.445906     1.130923\n",
       "Protein_1 3     I  0.159584     0.000000\n",
       "          15    V  0.168248     1.003839\n",
       "          22    L  0.589015     1.141511\n",
       "          24    P  0.523189     0.928291\n",
       "          25    D  0.804995     1.003516\n",
       "          27    V  1.404717     0.989936\n",
       "          29    M  1.662398     1.064119\n",
       "          32    L  1.071719     0.000000"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "comp = lambda x,y: x > y\n",
    "expressions=[(\"pcm\",comp,0)]\n",
    "\n",
    "merged_result.filter_result(expressions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Chapter 2: TAP prediction </h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fred2 offers only limited prediction methods for TAP prediction, due to lack of publicly available methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "svmtap 1.0\n",
      "smmtap 1.0\n",
      "doytchinova 1.0\n"
     ]
    }
   ],
   "source": [
    "for name,version in TAPPredictorFactory.available_methods().iteritems():\n",
    "    print name, \",\".join(version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For TAP prediction, we first have to generate peptides. Lets take the proteins we already imported and generate 9mers, as these two methods only support 9mer peptides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>svmtap</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(A, A, V, D, P, A, D, R, C)</th>\n",
       "      <td>-47.782823</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, D, R, C, K, E, V, Q, Q)</th>\n",
       "      <td>-38.337527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, V, D, P, A, D, R, C, K)</th>\n",
       "      <td>-40.900773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, G, K, A, A, V, D, P, A)</th>\n",
       "      <td>-31.988850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, K, E, V, Q, Q, I, R, D)</th>\n",
       "      <td>-38.903919</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                svmtap\n",
       "(A, A, V, D, P, A, D, R, C) -47.782823\n",
       "(A, D, R, C, K, E, V, Q, Q) -38.337527\n",
       "(A, V, D, P, A, D, R, C, K) -40.900773\n",
       "(C, G, K, A, A, V, D, P, A) -31.988850\n",
       "(C, K, E, V, Q, Q, I, R, D) -38.903919"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pep = list(generate_peptides_from_proteins(proteins,9))\n",
    "tap_result = TAPPredictorFactory(\"SVMTAP\").predict(pep)\n",
    "tap_result.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we can do all rudimentary operations on the result object as with the cleavage result objects, including merging and filtering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>svmtap</th>\n",
       "      <th>doytchinova</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(A, A, V, D, P, A, D, R, C)</th>\n",
       "      <td>-47.782823</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, D, R, C, K, E, V, Q, Q)</th>\n",
       "      <td>-38.337527</td>\n",
       "      <td>-0.466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, V, D, P, A, D, R, C, K)</th>\n",
       "      <td>-40.900773</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, G, K, A, A, V, D, P, A)</th>\n",
       "      <td>-31.988850</td>\n",
       "      <td>-2.070</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(C, K, E, V, Q, Q, I, R, D)</th>\n",
       "      <td>-38.903919</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                svmtap  doytchinova\n",
       "(A, A, V, D, P, A, D, R, C) -47.782823          NaN\n",
       "(A, D, R, C, K, E, V, Q, Q) -38.337527       -0.466\n",
       "(A, V, D, P, A, D, R, C, K) -40.900773          NaN\n",
       "(C, G, K, A, A, V, D, P, A) -31.988850       -2.070\n",
       "(C, K, E, V, Q, Q, I, R, D) -38.903919          NaN"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tap_result2 = TAPPredictorFactory(\"Doytchinova\").predict(pep[:15])\n",
    "tap_result.merge_results(tap_result2).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>svmtap</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(D, P, A, D, R, C, K, E, V)</th>\n",
       "      <td>-25.560005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(F, S, S, P, C, G, K, A, A)</th>\n",
       "      <td>-19.259415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(M, R, F, F, S, S, P, C, G)</th>\n",
       "      <td>-26.315791</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(S, P, C, G, K, A, A, V, D)</th>\n",
       "      <td>-24.671547</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(V, D, P, A, D, R, C, K, E)</th>\n",
       "      <td>-14.140251</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(V, P, D, H, V, N, M, S, E)</th>\n",
       "      <td>-23.865406</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                svmtap\n",
       "(D, P, A, D, R, C, K, E, V) -25.560005\n",
       "(F, S, S, P, C, G, K, A, A) -19.259415\n",
       "(M, R, F, F, S, S, P, C, G) -26.315791\n",
       "(S, P, C, G, K, A, A, V, D) -24.671547\n",
       "(V, D, P, A, D, R, C, K, E) -14.140251\n",
       "(V, P, D, H, V, N, M, S, E) -23.865406"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from operator import ge\n",
    "\n",
    "tap_result.filter_result([(\"svmtap\",ge, -30)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Chapter 3: Consensus prediction for natural ligand prediction </h2>\n",
    "<br/>\n",
    "Proteasomal cleavage, TAP prediction, as well as HLA binding prediction have been combined to increase the specificity of predicting natural processed HLA ligands. One example is `WAPP` <a href=\"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2279325/\">(Dönnes, et al.(2005). Protein Science 14(8))</a>, which uses proteasomal cleavage and TAP prediction methods to filter for possibly processed peptides.\n",
    "\n",
    "The same approach can be implemented with Fred2. We start again with the two protein sequences and exemplify the workflows for CleavageFragmentPrediction and CleavageSitePrediction methods.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4> Antigen processing prediction with CleavageFragment prediction </h4>\n",
    "<br/>\n",
    "We will use `PSSMGinodi`, `SVMTAP`, and `UniTope` for prediction for HLA-A*02:01. For `PSSMGionid` and `SVMTAP` we use a threshold of -15 and -30 respectively.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Number of peptides:  50\n",
      "Number of peptides after proteasomal cleavage:  31\n",
      "Number of peptides after TAP transport:  6\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>A*02:01</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Seq</th>\n",
       "      <th>Method</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(D, P, A, D, R, C, K, E, V)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-1.771554</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(F, S, S, P, C, G, K, A, A)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-1.221980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(M, R, F, F, S, S, P, C, G)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-2.002212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(S, P, C, G, K, A, A, V, D)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-5.042082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(V, D, P, A, D, R, C, K, E)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-4.721389</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(V, P, D, H, V, N, M, S, E)</th>\n",
       "      <th>unitope</th>\n",
       "      <td>-4.195872</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      A*02:01\n",
       "Seq                         Method           \n",
       "(D, P, A, D, R, C, K, E, V) unitope -1.771554\n",
       "(F, S, S, P, C, G, K, A, A) unitope -1.221980\n",
       "(M, R, F, F, S, S, P, C, G) unitope -2.002212\n",
       "(S, P, C, G, K, A, A, V, D) unitope -5.042082\n",
       "(V, D, P, A, D, R, C, K, E) unitope -4.721389\n",
       "(V, P, D, H, V, N, M, S, E) unitope -4.195872"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from operator import ge\n",
    "from Fred2.Core import Allele\n",
    "from Fred2.EpitopePrediction import EpitopePredictorFactory\n",
    "\n",
    "allele = Allele(\"HLA-A*02:01\")\n",
    "\n",
    "pep = list(generate_peptides_from_proteins(proteins,11))\n",
    "print \"Number of peptides: \", len(pep)\n",
    "\n",
    "#cleavage prediction and filtering\n",
    "df_cl = CleavageFragmentPredictorFactory(\"Ginodi\").predict(pep).filter_result((\"ginodi\",ge,-15))\n",
    "print \"Number of peptides after proteasomal cleavage: \", len(df_cl)\n",
    "\n",
    "#tap prediction and filtering\n",
    "df_tap = TAPPredictorFactory(\"SVMTAP\").predict(df_cl.index).filter_result((\"svmtap\",ge,-30))\n",
    "print \"Number of peptides after TAP transport: \", len(df_tap)\n",
    "\n",
    "#epitope prediction and filtering\n",
    "df_epi = EpitopePredictorFactory(\"UniTope\").predict(df_tap.index,alleles=allele)\n",
    "df_epi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on this analysis, there are no natural ligands predicted for the two test proteins."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4> Antigen processing prediction with CleavageSite prediction </h4>\n",
    "<br/>\n",
    "We will use `PCM`, `SVMTAP`, and `SVMHC` for prediction for HLA-A*02:01 like in the original work of WAPP. For `PCM` and `SVMTAP` we use a threshold of −4.8 and −27 respectively.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of peptides after proteasomal cleavage:  61\n",
      "Number of peptides after TAP transport:  5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>A*02:01</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Seq</th>\n",
       "      <th>Method</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(D, P, A, D, R, C, K, E, V)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-0.552764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(V, P, D, H, V, N, M, S, E)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-0.992899</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                     A*02:01\n",
       "Seq                         Method          \n",
       "(D, P, A, D, R, C, K, E, V) svmhc  -0.552764\n",
       "(V, P, D, H, V, N, M, S, E) svmhc  -0.992899"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from operator import ge\n",
    "from Fred2.Core import Allele,Protein\n",
    "from Fred2.EpitopePrediction import EpitopePredictorFactory\n",
    "\n",
    "allele = Allele(\"HLA-A*02:01\")\n",
    "\n",
    "\n",
    "\n",
    "#cleavage prediction and filtering\n",
    "df_cl = CleavageSitePredictorFactory(\"PCM\").predict(proteins).filter_result((\"pcm\",ge,-4.8))\n",
    "\n",
    "print \"Number of peptides after proteasomal cleavage: \", len(df_cl)\n",
    "\n",
    "#since we only predicted possible cleavage site, we now have to generate all possible peptides\n",
    "#peptides\n",
    "pep_dic = {}\n",
    "for p in proteins:\n",
    "    for i in df_cl.loc[(p.transcript_id,slice(None)),:].index.labels[1]:\n",
    "        if i-8>=0:\n",
    "            seq = str(p[i-8:i+1])\n",
    "            pep_dic.setdefault(seq, []).append(p)\n",
    "peps = [Peptide(seq, protein_pos={pp:[0] for pp in p}) for seq, p in pep_dic.iteritems()]\n",
    "\n",
    "\n",
    "#tap prediction and filtering\n",
    "df_tap = TAPPredictorFactory(\"SVMTAP\").predict(peps).filter_result((\"svmtap\",ge,-27))\n",
    "print \"Number of peptides after TAP transport: \", len(df_tap)\n",
    "\n",
    "#epitope prediction and filtering\n",
    "df_epi = EpitopePredictorFactory(\"SVMHC\").predict(df_tap.index,alleles=allele).filter_result((\"svmhc\",ge,-1.0))\n",
    "df_epi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here as well, we did not predict any natural ligands."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
